********************************************************************************
********************************************************************************
***                              Power Analysis                              ***
********************************************************************************
********************************************************************************
clear all

set matsize 800

*Create a matrix to store effect sizes
mat estimates = J(800,2,.)
mat pow = J(31,5,.)

*Number of subjects
//DKK 20.000, 11 minutes, 10.7 USD = DKK 74 minimum payment pr. hour, USD .25 = DKK 2 expected bonus payment
di 20000 / ((74+2) * 11/60) * 2/3

*Set number of observations
set obs 935
scalar subjects=935

*Create observations and randomly assign into treatment
gen id=_n
gen treatment=0
replace treatment=1 if id>subjects/2

**Loop over effect sizes
local countrow = 1
quietly forvalues i=4(.1)7 {
*Drop variables for each effect size (if there are any)
capture drop reg_pval MWU_pval significant_regtwo significant_MWUtwo significant_regone significant_MWUone

*Set treatment effect
scalar teffect = `i'

*Simulate 800 times
quietly forvalues j=1(1)800 {
*Drop values from previous loop (if there are any)
capture drop ftingroup ftoutgroup ftdiff

*Generate ft (incl. treatment effect), censor at 0 and 100, and examine difference
gen ftingroup = rnormal(62, 24)
replace ftingroup = 0 if ftingroup<0
replace ftingroup = 100 if ftingroup>100

gen ftoutgroup = rnormal(40, 26)
replace ftoutgroup = 0 if ftoutgroup<0
replace ftoutgroup = ftoutgroup - teffect if treatment==1
replace ftoutgroup = 100 if ftoutgroup>100

gen ftdiff = ftingroup - ftoutgroup

*Saving p-value from regression
reg ftdiff treatment, robust
matrix estimates[`j',1]=r(table)[4,1]

*Saving p-value from Mann-Whitney U-test
ranksum ftdiff, by(treatment)
matrix estimates[`j',2]=r(p)
}

*Making variables of the matrix
svmat estimates, names(temp)

*Rename variables
rename temp1 reg_pval
rename temp2 MWU_pval

*Save data set for later use
*save $input/powerbetween.dta, replace

*Compare power of tests
gen significant_regtwo = 0 if reg_pval!=.
replace significant_regtwo = 1 if reg_pval<0.05
gen significant_MWUtwo = 0 if MWU_pval!=.
replace significant_MWUtwo = 1 if MWU_pval<0.05

gen significant_regone = 0 if reg_pval!=.
replace significant_regone = 1 if reg_pval<0.1
gen significant_MWUone = 0 if MWU_pval!=.
replace significant_MWUone = 1 if MWU_pval<0.1

*Save power in a matrix
matrix pow[`countrow',1]=`i'

ci means significant_regtwo
matrix pow[`countrow',2]=r(mean)

ci means significant_MWUtwo
matrix pow[`countrow',3]=r(mean)

ci means significant_regone
matrix pow[`countrow',4]=r(mean)

ci means significant_MWUone
matrix pow[`countrow',5]=r(mean)

local countrow = `countrow'+1
*Showing simulation number
noisily display `i'
}

*Making variables of the power matrix
svmat pow, names(temp)

*Rename variables
rename temp1 effectsize
rename temp2 reg_powertwo
rename temp3 MWU_powertwo
rename temp4 reg_powerone
rename temp5 MWU_powerone

*Calculate minimum detectable effect size
esizei 475 22 35 475 16.6 35
esizei 475 22 35 475 16 35
